单细胞转录组高级分析五：GSEA与GSVA分析

上期专题我们介绍了单细胞转录组数据的基础分析，然而那些分析只是揭开了组织异质性的面纱，还有更多的生命奥秘隐藏在数据中等待我们发掘。本专题将介绍一些单细胞转录组的高级分析内容：多样本批次校正、SCENIC转录调控网络分析、CellPhoneDB细胞通讯分析、inferCNV肿瘤拷贝数变异分析、GSEA/GSVA分析和其他高级分析内容。不足之处请大家批评指正，欢迎添加Kinesin微信交流探讨！（扫描文末二维码）

往期相关文章

单细胞转录组基础分析专题

GSEA与GSVA简介
基因集的概念

GSEA全称Gene Set Enrichment Analysis，GSVA全称Gene Set Variation Analysis，它们都是基于基因集开展的分析，因此我们先要了解基因集的定义。基因集顾名思义就是一些基因的集合，任何一些基因放在一起都可以叫做基因集，但是我们用来分析的基因集要求有一定的生物学意义。最常见的基因集数据库如GO与KEGG，它们一个按照基因本体论将基因分门别类，一个按照代谢通路将相关基因集合在一起。除此之外我们还可以按转录因子调控网络、共表达网络、定义生物状态的marker基因列表等把基因集合成有一定生物学意义的基因集。

MSigDB基因集数据库

GSEA是由Broad研究所开发的一种富集方法，他们在提出该方法的同时还提供了一个基因集数据库——MSigdb。它从位置，功能，代谢途径，靶标结合等多种角度出发，构建出了许多的基因集合，Broad研究所将他们构建的基因集合保存在MSigDB，官网地址如下：

http://software.broadinstitute.org/gsea/msigdb/index.jsp

在MSigDB中，将所有的基因集划分为以下8大类别：

1. H：hallmark gene sets

特征基因集，由定义生物状态和进程的marker基因组成。

2. C1：positional gene sets

位置基因集，包含人类每条染色体上的不同cytoband区域对应的基因集合。

3. C2：curated gene sets

代谢通路基因集，包含KEGG, Reactome, BioCarta数据库，以及文献和专家支持的基因集信息。

4. C3：motif gene sets

靶基因集，包含了miRNA靶基因集和转录因子调控基因集两大类。

5. C4：computational gene sets

计算基因集，计算机软件预测出来的基因集，主要是和癌症相关的基因。

6. C5：GO gene sets

基因本体基因集，包含了Gene Ontology对应的基因集合。

7. C6：oncogenic signatures

癌症扰动基因集，来源于药物处理肿瘤后基因差异表达数据，包含已知条件处理后基因表达量发生变化的基因。

8. C7：immunologic signatures

免疫基因集，包含了免疫系统功能相关的基因集合。

GSEA的分析原理

常规GO/KEGG富集分析需要设定阈值过滤差异基因，阈值太宽富集的结果太多，阈值太严又可能会遗漏一些关键结果。GO/KEGG富集的结果通常还很宽泛，并不能很好地解释生物学现象。有鉴于此，Broad研究所开发了基因集富集分析(GSEA)方法。GSEA使用无监督算法，不用过滤任何基因，配合MSigDB数据库使用，更容易找到解释生物学现象的基因集。其原理如下：

GSEA分析要先将样本做组间对比分析，GSEA自带9种分析方法，分为基因表达值差异分析和相关性分析两大类。对于Case/Contral的实验设计，差异分析方法更为常用，这其中又以默认的信噪比和大家熟悉的差异倍数用的最多。对比分析之后要按结果将基因排序，以差异倍数方法为例，把所有基因按差异倍数(FC)的值降序排列以供后续分析。上图小人脚下的小方块代表排序好的差异基因列表，蓝色之外的其他色块代表属于某个基因集的基因，如黄色属于基因集A，绿色属于基因集B。最下面高低不等的竖条代表与基因列表对应的FC值，红色上调、蓝色下调、黄色没有变化。基因集的富集分析需要经历三步：

1. 基因集A富集分析时，小人从基因列表的左端走到右端，每经过一个蓝色基因扣分，每遇到一个黄色基因加分，扣分时与FC无关，加分时考虑FC的权重。基因集A最终的富集分数(ES)是小人曾经得过的最高/低分，实际公式比这复杂，但基本理念如此。

2. 采用置换检验计算基因集A的显著性，即p值。

3. 基因集A富集分析完成后，按上述同样的方法完成基因集B、C直至所有输入基因集的分析。所有需要富集分析的基因集都计算ES和p值之后，将ES转换为标准富集分数(NES)，并计算校正后p值。
   

听完我的解释之后再看官方的解释可能更容易理解：

A GSEA overview illustrating the method. (A) An expression dataset sorted by correlation with phenotype, the corresponding heat map, and the ‘‘gene tags,’’ i.e., location of genes from a set S within the sorted list. (B) Plot of the running sum for S in the dataset, including the location of the maximum enrichment score (ES) and the leading-edge subset.

GSVA的分析原理

GSEA虽然是一种强大的富集分析工具，但是它的应用场景通常局限于Case/Control的实验设计。对于表型（分组）复杂的大样本量研究，例如TCGA和单细胞图谱这样的项目，分析起来就非常困难。因此，Broad研究所在GSEA发布8年之后，又开发了GSVA算法来拓展基因集分析的应用。GSVA不需要预先进行样本之间的差异分析，它依据表达矩阵就可以计算每个样本中特定基因集的变异分数。简单的说，输入以基因为行的表达矩阵和基因集数据库给GSVA，它就输出以基因集名称为行的变异分数矩阵，如下图所示：

左侧输入基因表达矩阵和基因集数据库，中间是GSVA算法原理，右侧是输出的基因集变异分数矩阵。基因集变异分数可以理解为基因集内所有基因的综合表达值。

GSEA的安装与使用

GSEA官方提供很多版本，但是不提倡R语言版，下载网址如下：

https://www.gsea-msigdb.org/gsea/downloads.jsp

准备分析文件

单细胞GSEA分析的难点是输入文件的准备，大家在R语言中运行下列代码，就可以得到符合GSEA要求的特定格式文件。

##R语言准备gsea输入文件
library(Seurat)
library(tidyverse)
library(GSEA)
dir.create("GSEA")   
dir.create("GSEA/input")
dir.create("GSEA/output")
scRNA <- readRDS("scRNA.rds")
#提取HNC01TIL和Tonsil2的T细胞做富集分析
tmp <- subset(tmp, subset = (tmp$orig.ident=='HNC01TIL'|tmp$orig.ident=='Tonsil2')&tmp$celltype_Monaco=='T cells')
sub.cells <- rownames(tmp)
scRNAsub <- subset(scRNA, cells=sub.cells)
expr <- GetAssayData(scRNAsub, slot = 'counts')
expr <- data.frame(NAME=rownames(expr), Description=rep('na', nrow(expr)), expr, stringsAsFactors=F)
write('#1.2', "GSEA/input/expr.gct", ncolumns=1)
write(c(nrow(expr),(ncol(expr)-2)), "GSEA/input/expr.gct", ncolumns=2, append=T, sep='\t')
write.table(expr, "GSEA/input/expr.gct", row.names=F, sep='\t', append=T, quote=F)
line.1 <- c((ncol(expr)-2), 2, 1)
tmp <- table([e-mail]$orig.ident)
line.2 <- c("#", names(tmp))
line.3 <- c(rep(names(tmp)[1],tmp[1]), rep(names(tmp)[2],tmp[2]))
write(line.1, 'GSEA/input/group.cls', ncolumns=length(line.1), append=T, sep='\t')
write(line.2, 'GSEA/input/group.cls', ncolumns=length(line.2), append=T, sep='\t')
write(line.3, 'GSEA/input/group.cls', ncolumns=length(line.3), append=T, sep='\t')

准备好的文件格式：

表达数据gct文件内容节选

分组cls文件内容节选

得到上面的gct和cls文件，再加上基因集数据库gmt文件，就可以进行GSEA分析了。建议在java版本上运行，官方的R包太简陋了，java版上的好多功能都不支持。

R语言版代码与结果展示

##安装官方R包GSEA
library(devtools)
install_github("GSEA-MSigDB/GSEA_R")
#gs.db参数对应的是基因集数据文件，这里只对kegg收录的基因集做了分析
GSEA::GSEA("GSEA/input/expr.gct", "GSEA/input/group.cls", 
           gs.db="~/database/GSEA/c2.cp.kegg.v7.1.symbols.gmt", output.directory="GSEA/output/")

分析得到的结果很多，分为结果总览，分组概况，以及各个基因集的富集情况详细。这里展示NES分数最高的通路的富集结果：

上图对应数据

java版分析代码

提前安装好java11，下载图中所示GSEA版本，解压即可

#以下命令在linux环境中运行，windows环境用gsea-cli.bat代替gsea-cli.sh
bash gsea-cli.sh GSEA -res ~/project/2020/2007_10xDemo2/GSEA/input/expr.gct \
                      -cls ~/project/2020/2007_10xDemo2/GSEA/input/group.cls#HNC01TIL_versus_Tonsil2 \
                      -gmx ~/project/2020/2007_10xDemo2/GSEA/input/c2.cp.kegg.v7.1.symbols.gmt \
                      -collapse false \
                      -out ~/project/2020/2007_10xDemo2/GSEA/out/

图形界面版GSEA操作

表达数据和分组文件准备好之后，图形界面版操作没有难度，简单介绍一下：

这是开始界面，先加载数据

点击箭头所指按钮之后，就会弹出对话框让你选择文件，分次输入我们之前准备好的表达数据文件(*.gct)和分组文件(*.cls)

选好文件之后，文件目录会出现在左右两边箭头所指的地方。如果右边没有，双击左边即可。我这儿左边有三个文件是因为加载了本地保存的基因集数据库文件。

文件加载完毕之后，点击箭头所指运行会出现右边的对话框，注意一下红框标注的选项。Collapse/Remap项选择“No_Collapse”（Collapse适用芯片数据），Permutation项选择“phenotype”（更适合单细胞数据），Chip项可不选或选择有“Human_Symbol_with_”关键字的项目（选择之后结果中基因有简短注释）。

Basic fields点开之后选择保存结果的目录，点击“Run”就开始运行了。java版（图形界面的也是java版）运行之后结果与R语言版一致，不过分析报告是网页版，点击超链接就可以查看细节，使用起来非常方便。

图形也更美观

GSVA的安装与使用

安装GSVA

#从bioconductor安装
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("GSVA")
#或者从github安装
library(devtools)
install_github("rcastelo/GSVA")

准备输入数据

GSVA要求输入表达矩阵和基因集列表。表达矩阵从seurat对象导入即可，以列表的形式提供基因集数据有一定难度，我编写了一个函数解决此问题。

library(Seurat)
library(GSVA)
library(tidyverse)
##创建gmt文件转list函数
gmt2list <- function(gmtfile){
 sets <- as.list(read_lines(gmtfile))
 for(i in 1:length(sets)){
    tmp = str_split(sets[[i]], '\t')
  n = length(tmp[[1]])
  names(sets)[i] = tmp[[1]][1]
  sets[[i]] = tmp[[1]][3:n]
  rm(tmp, n)
 }
 return(sets)
}
#读取基因集数据库
s.sets = gmt2list("GSEA/input/c2.cp.kegg.v7.1.symbols.gmt")
#读取表达矩阵
scRNA <- readRDS("scRNA.rds")
# 随机提取1000个细胞演示GSVA，非常规操作
# tmp <- sample(colnames(scRNA),1000) %>% sort()
# scRNA <- scRNA[,tmp]
expr <- as.matrix(scRNA@assays$RNA@counts)
meta <- [e-mail][,c("seurat_clusters", "celltype_Monaco")]
es.matrix = gsva(expr, s.sets, kcdf="Poisson")
write.table(es.matrix, 'GSVA/gsva.xls', row.names=T, col.names=NA, sep='\t')

结果可视化

GSVA的分析结果相当于表达矩阵，用热图来可视化比较常用，为了体现单细胞的特点我也用FeaturePlot展示了几个基因集。

library(pheatmap)
library(patchwork)
#绘制热图
pheatmap(es.matrix, show_rownames=1, show_colnames=0, annotation_col=meta,
         fontsize_row=5, filename='GSVA/gsva_demo.png', width=15, height=12)
#挑选感兴趣的基因集绘制featureplot
es <- data.frame(t(es.matrix),stringsAsFactors=F)
scRNA <- AddMetaData(scRNA, es)
p1 <- FeaturePlot(scRNA, features = "KEGG_PRIMARY_BILE_ACID_BIOSYNTHESIS", reduction = 'tsne')
p2 <- FeaturePlot(scRNA, features = "KEGG_ETHER_LIPID_METABOLISM", reduction = 'tsne')
p3 <- FeaturePlot(scRNA, features = "KEGG_RIBOSOME", reduction = 'tsne')
p4 <- FeaturePlot(scRNA, features = "KEGG_ASTHMA", reduction = 'tsne')
plotc = (p1|p2)/(p3|p4)
ggsave('GSVA/featureplot_demo.png', plotc, width = 10, height = 8)

随意挑了几个基因集展示gsva富集分数

GSVA分析给我们的不是阶段性结果，而是观察数据的一个新的视角，让我们以基因集为单位研究数据。基于基因集的富集矩阵，我们可以做差异分析、表型相关性分析、生存分析，也可以基于基因集的功能来鉴定细胞类型。

References

https://www.gsea-msigdb.org/gsea/doc/GSEAUserGuideFrame.html

https://www.gsea-msigdb.org/gsea/msigdb/collection_details.jsp

Laurila E , Mootha V K , Lindgren C M , et al. PGC-1α-responsive genes involved in oxidative phosphorylation are coordinately downregulated in human diabetes[J]. 2003, 34(3):267-273.

Subramanian, Tamayo, et al. Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles. 2005, PNAS 102, 15545-15550.

Hnzelmann S , Castelo R , Guinney J . GSVA: gene set variation analysis for microarray and RNA-Seq data[J]. Bmc Bioinformatics, 2013, 14(1):7-7.

获取帮助

本教程的目的在于把单细胞分析流程串起来，适合有一定R语言基础的朋友参考。分析原理和代码我没有详细解释，官网有详细的教程和权威的解释，翻译好的中文教程也有多个版本，有兴趣的可以搜索一下。如果您学习本教程有一定困难，可以分享此篇文章到朋友圈，截图微信发给Kinesin（文末二维码），我会抽时间组群直播讲解单细胞数据分析的全过程。本专题所用的软件、数据、代码脚本和中间结果，我打包放在了百度云上，需要的朋友添加Kinesin微信索取。